Intersections module

Set up plotly


In [1]:
import plotly.offline as py
import plotly.graph_objs as go
py.init_notebook_mode(connected=True)


Set up imports


In [2]:
import pointslinesplanes as plp
import numpy as np

In [3]:
def intersect(obj1, obj2):
    """Intersection between two objects that can be Point, Line or Plane.
    Returns Point, Line or Plane or throws ArithmeticError when no intersection can be found.
    @author Nick Metelski
    @since 26.07.17"""
    #plane-plane
    if isinstance(obj1,plp.Plane) and isinstance(obj2,plp.Plane):
        return _plane_plane_intersect(obj1,obj2)
    
    #plane-line
    elif isinstance(obj1,plp.Plane) and isinstance(obj2,plp.Line):
        return _plane_line_intersect(obj1,obj2)
    elif isinstance(obj1,plp.Line) and isinstance(obj2,plp.Plane):
        return _plane_line_intersect(obj2,obj1)
    
    #plane-point
    elif isinstance(obj1,plp.Plane) and isinstance(obj2,plp.Point):
        return _plane_point_intersect(obj1,obj2)
    elif isinstance(obj1,plp.Point) and isinstance(obj2,plp.Plane):
        return _plane_point_intersect(obj2,obj1)
    
    #line-line
    elif isinstance(obj1,plp.Line) and isinstance(obj2,plp.Line):
        return _line_line_intersect(obj1,obj2)
    
    #line-point
    elif isinstance(obj1,plp.Line) and isinstance(obj2,plp.Point):
        return _line_point_intersect(obj1,obj2)
    elif isinstance(obj1,plp.Point) and isinstance(obj2,plp.Line):
        return _line_point_intersect(obj2,obj1)
    
    #point-point
    elif isinstance(obj1,plp.Point) and isinstance(obj2,plp.Point):
        return _point_point_intersect(obj1,obj2)
    
    #wrong params
    else:
        raise TypeError("Invalid parameter types - please pass intersections.Point, intersections.Line or intersections.Plane.")

def _plane_plane_intersect(p1,p2):
    """Private function; do not use. Use intersect(obj1,obj2) instead.
    plane-plane intersection submethod. Raises exceptions or returns Line.
    @author Nick Metelski
    @since 26.07.17"""
    #TODO make it return plane is planes are identical
    #cross-product
    cross = np.cross(p1.normal,p2.normal)
    if np.all(cross==0):
        #planes are parallel or overlap
        #if they overlap then the offset of one plane lies on the other
        try:
            intersection = intersect(Point(p1.offset), p2)
            return Plane(p1.normal,p1.offset) #return a new plane same a the overlapping ones
        except ArithmeticError:
            raise ArithmeticError("Planes are parallel, do not overlap.")
    else:
        #sample point: x=0
        #NOTE: linalg.solve can raise np.linalg.LinAlgError exception when x=0 results in singular matrix
        #we have to check some other cases
        #premise: the line has to intersect at least one of the planes: xy, yz or xz.
        mat = [[p1.normal[1],p1.normal[2]],[p2.normal[1], p2.normal[2]]]
        axis = 0 #keep track of which axis we zero out; 0=x, 1=y, 2=z
        if np.linalg.matrix_rank(mat) == 1:
            mat = [[p1.normal[0],p1.normal[2]],[p2.normal[0], p2.normal[2]]]
            axis = 1
            if np.linalg.matrix_rank(mat) == 1:
                mat = [[p1.normal[0],p1.normal[1]],[p1.normal[0], p2.normal[1]]]
                axis = 2  
        rhs = [np.dot(p1.normal,p1.offset),np.dot(p2.normal,p2.offset)]
        sol = np.linalg.solve(mat,rhs)
        if axis == 0:
            return plp.Line(cross, [0,sol[0],sol[1]])
        if axis == 1:
            return plp.Line(cross, [sol[0],0,sol[1]])
        if axis == 2:
            return plp.Line(cross, [sol[0],sol[1],0])

def _plane_line_intersect(plane,line):
    """Private function; do not use. Use intersect(obj1,obj2) instead.
    Intersection of a line an a plane. Raises exceptions or returns Point.
    @author Nick Metelski
    @since 26.07.17"""
    check = np.dot(plane.normal,line.vec)
    if np.allclose(check,[0.0]):
        #plane and line are parallel or overlap
        raise ArithmeticError("Plane and line are parallel.")
    else:
        #there is an explicit formula for parameter of the line for which intersection is met:
        #$$t = \frac{\vec{n} \cdot (\vec{d_p}-\vec{d_v})}{\vec{n} \cdot \vec{v}}$$,
        #where n is normal to plane, v is line's direction, rp is plane offset, rv is vector offset
        t = np.dot(plane.normal,(plane.offset-line.offset))/np.dot(plane.normal,line.vec)
        return plp.Point(line.offset+t*line.vec)
    
def _plane_point_intersect(plane,point):
    """Private function; do not use. Use intersect(obj1,obj2) instead.
    Intersection of a plane and point. Raises exceptions or returns Point.
    @author Nick Metelski
    @since 27.07.17"""
    check = np.dot(plane.normal, point.pos - plane.offset)
    if check == 0:
        return plp.Point(point.pos)
    else:
        raise ArithmeticError("The point does not line on the plane.")

def _line_line_intersect(line1,line2):
    """Private function; do not use. Use intersect(obj1,obj2) instead.
    Intersection of two lines. Raises exceptions or returns Point or Line.
    @author Nick Metelski
    @since 27.07.17"""
    #raise NotImplementedError("Line-line intersections not implemented.")
    cross = np.cross(line1.vec,line2.vec) #is normalized as vectors are normalized
    if np.allclose(cross,[0.0,0.0,0.0]):
        #lines parallel. are they identical?
        try:
            intersect(line2, plp.Point(line1.offset)) #this will raise an exception if they are not identical
            return plp.Line(line1.vec,line1.offset)
        except ArithmeticError:
            raise ArithmeticError("Lines are parallel, and are not identical.")
    #lines not parallel
    sample = line2.offset - line1.offset #sample vector between two points on the line
    dist = np.abs(np.dot(sample,cross)) #distance between the lines
    if np.allclose([dist],[0.0]):
        v1 = line1.vec
        v2 = line2.vec
        d1 = line1.offset
        d2 = line2.offset
        n = plp.utils.normalize(np.cross(v2,v1))
        s = d2 - d1
        l = np.linalg.norm(np.cross(s,v1))
        l_vec = l * np.cross(v1,n)
        param = l**2/np.dot(l_vec,v2)

        return plp.Point(line2.offset+param*line2.vec)
    else:
        raise ArithmeticError("Lines are skew, no intersection. Distance: " + str(dist))
    
def _line_point_intersect(line,point):
    """Private function; do not use. Use intersect(obj1,obj2) instead.
    Intersection of a line and point. Intersection is a Point.
    Raises exceptions or returns Point.
    @author Nick Metelski
    @since 27.07.17"""
    #premise: difference in offsets must be parallel to direction vector
    offset_diff = point.pos - line.offset
    cross = np.cross(line.vec,offset_diff) # this is for checking if there exists an intersection
    if np.allclose(cross,[0.0,0.0,0.0]):
        return plp.Point(point.pos) #return a point coinciding with the original point passed as param
    else:
        raise ArithmeticError("Point does not lie on the line.")
        
def _point_point_intersect(point1, point2):
    """Private function; do not use. Use intersect(obj1,obj2) instead.
    Intersection of a two points. Intersection is the point if the point are the same position. 
    Raises exceptions or returns Point.
    @author Nick Metelski
    @since 27.07.17"""
    if np.allclose(point1.pos,point2.pos):
        #just return one of the points when they coincide
        return plp.Point(point1.pos)
    else:
        raise ArithmeticError("Points do not coincide.")

Tests


In [4]:
plane1 = plp.Plane([1,1,1],[1,0,0])
plane2 = plp.Plane([-1.0,1,1],[0.5,0,0.5])
plane3 = plp.Plane([0.2,0.1,0.3],[0.2,-0.3,0.5])
line1 = intersect(plane1,plane2)
line2 = intersect(plane2,plane3)
line3 = intersect(plane1,plane3)
point1 = intersect(line1,line2)

In [5]:
layout = dict(
    width=1000,height=1000,
    showlegend=False,
    font = dict(family="Verdana"),
    scene = dict(
        xaxis = dict(range=[-1, 1], autorange=False, zeroline=False),
        yaxis = dict(range=[-1, 1], autorange=False, zeroline=False),
        zaxis = dict(range=[-1, 1], autorange=False, zeroline=False),
        aspectmode = 'cube',
        camera = dict(center=dict(x=0,y=0,z=0),eye=dict(x=0,y=0,z=3))
    ),
    plot_bgcolor='rgb(255, 255, 255)'
)

In [6]:
objs = [plane1, plane2, plane3, line1, line2, line3, point1]
data = [o.goify() for o in objs]
fig = go.Figure(data=data,layout=layout)
plot = py.iplot(fig)



In [29]:
#Line(dir=Vector(0.447,0.894,0.000),off=Vector(0.120,-0.060,0.000)) and Line(dir=Vector(0.316,0.949,0.000),off=Vector(0.150,-0.050,0.000))
line2 = plp.Line([0.4472135954999579,0.8944271909999159,0],[0.120,-0.06,0])
line1 = plp.Line([0.31622776601683794,0.9486832980505138,0],[0.15,-0.05000000000000002,0])

v1 = line1.vec
v2 = line2.vec
print("v2",v2)
d1 = line1.offset
d2 = line2.offset
n = plp.utils.normalize(np.cross(v2,v1))
print(n)
s = d2 - d1
print(s)
l = np.linalg.norm(np.cross(s,v1))
print(l)
l_vec = l * np.cross(v1,n)
print("v1 x n", np.cross(v1,n))
print("l_vec ", l_vec)
print("lxl", l*l)
print("dot",np.dot(l_vec,v2))
param = (l*l)/np.dot(l_vec,v2)
print(param)

point = plp.Point(line2.offset+param*line2.vec)


v2 [ 0.4472136   0.89442719  0.        ]
[ 0.  0.  1.]
[-0.03 -0.01  0.  ]
0.0252982212813
v1 x n [ 0.9486833  -0.31622777  0.        ]
l_vec  [ 0.024 -0.008  0.   ]
lxl 0.00064
dot 0.003577708764
0.1788854382

In [28]:
objs = [line1,line2,point]
data = [o.goify() for o in objs]
fig = go.Figure(data=data,layout=layout)
plot = py.iplot(fig)
print(point.goify())


{'type': 'scatter3d', 'x': [0.20000000000000001], 'mode': 'markers', 'y': [0.10000000000000001], 'z': [0.0]}
$$\vec{r} = \vec{d_2} + \Bigg(\Big[ \big| (\vec{d_2} - \vec{d_1})\times \vec{v_1} \big| \bullet \big( \vec{v_1} \times ( \vec{v_2} \times \vec{v_1}) \big) \Big] \cdot \vec{v_2}\Bigg) \bullet \vec{v_2}$$

In [34]:
plane = plp.Plane([0.1,0.3,-0.2],[0.2,0.1,0.4])
line = plp.Line([0.1,0.3,0],[0.2,0.1,0])
pt = intersect(plane,line)

In [35]:
objs = [plane,line,pt]
data = [o.goify() for o in objs]
fig = go.Figure(data=data,layout=layout)
plot = py.iplot(fig)
print(point.goify())


{'type': 'scatter3d', 'x': [0.20000000000000001], 'mode': 'markers', 'y': [0.10000000000000001], 'z': [0.0]}

In [ ]: